% This script starts estimation of the Macro finance model using the SR approach
% By Martin M. Andreasen, August 2019
addpath(genpath(pwd))
close all
clear 
%clc
%% User settings
dataOption        = 1;          %1 for using ygap = UGAP,
                                %2 for using ygap = potential output in GDP
                                %3 for using ygap = Cooper and Priestley gap
startYear         = 1961;       %startYear for the estimation. Data available from 1961
modelNum          = 1;          %1 for QTSM without interest rate smoothing
                                %2 for QTSM wiht constant interest rate smoothing

% For inference
theta12OptimalOn  = 0;          %1 for setting LAMBDA in step 3 in an optimal manner. 0 for using the estimates from step 2
AVarTheta1On      = 0;          %1 for computing AVar for theta1 at step 1 and theta11 at step 3, else 0

% The optimizer:
optimizerStep1     = 0;         % optim = 0 for no optimization
optimizerStep3     = 0;         % optim = 1 for the CMAES
                                % optim = 2 for a gradient based optimizer (the LM version)
                                % optim = 3 for the Nelder-Mead.
                                % optim = 4, for first a gradient optimizer and the the Nelder-Mead
optim.numOptim    = 1;          % Number of optimizations - if optimizer = 3 or 4.
optim.sigma       = 0.10;       % step length for CMAES
optim.sigmaMax    = 0.20;       % max length for CMAES
optim.TolFun      = 1e-5;
optim.TolX        = 1e-5;
optim.MaxIter     = 4000;
optim.MaxFunEvals = 4000;

epsValue          = 1D-6;       % Stepsize for the standard errors of Theta1
AVarAdjThetaOn    = 0;          % 1 to adjust factor variance for Theta1, else 0

h0RegimeOn        = 0;          %1 for switching in h0 under P, else 0
hxRegimeOn        = 1;          %1 for switching in hx under P, else 0

%% Default setting
numCPUs           = feature('numcores'); % Number of CPUs
fullPHIon         = 1;          %1 for a full phi, else 0 for diagonal phi
nm                = 2;          %Number of macro factors
nn                = 0;          %Number of latent factors factors
loadDefaultSetting

% Loading the data:
loadData = getData(dataOption,startYear);
%% Starting values for the estimated parameters
theta1 = initializeTheta1(fullPHIon,nm,nn);

% Impose restrictions on phi from ReStud paper
theta1 = rmfield(theta1,'phi14');
theta1 = rmfield(theta1,'phi23');
theta1 = rmfield(theta1,'phi32');
theta1 = rmfield(theta1,'phi34');
theta1 = rmfield(theta1,'phi41');
theta1 = rmfield(theta1,'phi43');

if modelNum == 1  %
    theta1 = rmfield(theta1,'rhor');
    theta1ValuesStep1 = [   0.000007697957    0.088798954825    0.000001448479    0.006434745357    0.025271178113  ...
        0.003801472813    0.004258045841   -0.485207099637   -0.491153447384   -0.021787206525  ...
        -0.000928735436    0.002855108937   -0.024156823397    0.000103525349    0.000830331327  ...
        0.023919846878   -0.087674048434    0.049990615980    0.018135323127   -0.000302856709 ...
        0.039970173130    0.002512837718   -0.004243300828    0.002117672196    0.000003903177  ]';
    % Q = 0.069862979582060, done by CMAES
    % For step 3
    if h0RegimeOn == 0 && hxRegimeOn == 0
        theta1ValuesStep3 = [   -0.000054267465    0.062141749544    0.000001029203    0.020190705171    0.023033327077  ...
            0.002630766803    0.004823063396   -0.252093024059   -0.361085717407    0.015628416992  ...
            -0.000725748632    0.002047542392   -0.015202204828    0.004076399520    0.000672800837  ...
            0.001655995736    0.000039057451    0.001943359316   -0.006279208110   -0.001922424891  ...
            0.011801466352   -0.000719350910   -0.000207358753    0.000926326930    0.000712422918 ]';
        % Q =   0.073725835300225, done by CMAES
    elseif h0RegimeOn == 1 && hxRegimeOn == 0
        theta1ValuesStep3 = [ ]';
        % Q =
    elseif h0RegimeOn == 0 && hxRegimeOn == 1
        theta1ValuesStep3 = [   -0.000052529151    0.061893644107    0.000001015031    0.020127398256    0.023015775710  ...
            0.002649581929    0.004818868026   -0.249851076249   -0.362188831135    0.015830591983  ...
            -0.000727219081    0.002044829664   -0.015189659630    0.004084626447    0.000671770332  ...
            0.001634667604    0.000033630677    0.001939529654   -0.006542854015   -0.001840507312 ...
            0.011472572799   -0.000733829336   -0.000202007946    0.000909291558    0.000711139421 ]';
        % Q= 0.073727371238604, done 
    elseif h0RegimeOn == 1 && hxRegimeOn == 1
        % Change in intercepts are not significant
        theta1ValuesStep3 = [     -0.000105253887    0.414352586444    0.494286583897    0.301285280938    0.000001093386  ...
            0.073390765022    0.287349341969   -0.041571247513   -0.003714475854    0.389062513331  ...
            -0.088959177209    0.037238279387    0.085480399328    0.026462625781   -0.015214482809  ...
            0.001693778664    0.000059919995    0.001949913807   -0.006211423343   -0.001786126762 ....
            0.011687749617   -0.000721208403   -0.000199725004    0.000921912262    0.000711221986  ]';
        % Q = 0.096928521607404, done by CMAES
    end
elseif modelNum == 2
    
end
theta1           = values2struct(theta1ValuesStep1,fieldnames(theta1));
theta1StartStep3 = values2struct(theta1ValuesStep3,fieldnames(theta1));

% Getting bounds and Insigma
[upperBounds,lowerBounds,Insigma] = theta1BoundsInsigma(nm,nn);

% Defining the struct with the setup specification for step 1
constructSetupStep1

%% Computing the objective function theta1
tic
[Q,output,errorMes] = objectFunc(struc2values(theta1),setupStep1);
toc
%% Step 1 of the SR approach
% The optimization
optim.optimizer = optimizerStep1;
[QStep1,theta1Step1,outputStep1] = step1SR(theta1,setupStep1,optim);
theta1 = theta1Step1;

% The standard errors
AVarTheta1Step1  = getAVarTheta1andFactors(theta1Step1,setupStep1,epsValue);
%plotFactors(setupStep1,outputStep1,AVarTheta1Step1)
%% Step 2 of the SR approach
[theta2Step2,AVarTheta2Step2] = step2SR_threshold_macro(AVarTheta1Step1,outputStep1,bandwidthTstep2,setupStep1.econAct,h0RegimeOn,hxRegimeOn,numBootStep2);
%mprStep2 = getMPR(outputStep1,theta2Step2,AVarTheta2Step2,nm*2+nn,h0RegimeOn,hxRegimeOn);
%% Step 3 of the SR approach
namesTheta12 = cell((2*nm+nn)*(2*nm+nn+1)/2,1);
idx = 0;
for i=1:2*nm+nn
    for j=1:i
        idx = idx + 1;
        namesTheta12{idx} = ['sigma',num2str(i),num2str(j)];
    end
end
[theta12Step3,AVarTheta12Step3] = step3SR_theta12(namesTheta12,theta1Step1,AVarTheta1Step1,theta2Step2,AVarTheta2Step2,theta12OptimalOn);

% Estimation of theta11 - using Robust version
if exist('theta1StartStep3','var')
    theta1Start = theta1StartStep3;
    if AVarTheta1On  == 1
        disp('for standard errors')
        theta1Start = rmfield(theta1Start,'phi22');
        setupStep1.calibrateTheta1.phi22 = theta1StartStep3.phi22;
        setupStep1.selectTheta1          = fieldnames(theta1Start);       
        setupStep1.upperBoundsValues    = struc2values(upperBounds,setupStep1.selectTheta1);
        setupStep1.lowerBoundsValues    = struc2values(lowerBounds,setupStep1.selectTheta1);
        setupStep1.InsigmaValues        = struc2values(Insigma,setupStep1.selectTheta1);
     end
else
    theta1Start = theta1Step1;
end
optim.optimizer = optimizerStep3;
[QStep3,theta11Step3,outputStep3,setupStep3] = step3SR_theta11(theta12Step3.Robust,namesTheta12,theta1Start,...
    setupStep1,upperBounds,lowerBounds,Insigma,optim);

% Standard errors of theta11: The asymptotic standard errors 
AVarTheta11Step3 = getAVarTheta11Step3(namesTheta12,theta11Step3,AVarTheta1Step1,AVarTheta2Step2,AVarTheta12Step3,...
    setupStep1,setupStep3,epsValue,theta12OptimalOn);

% Re-estimating P dynamics
[theta2Step3,AVarTheta2Step3] =step2SR_threshold_macro(AVarTheta11Step3,outputStep3,bandwidthTstep2,setupStep3.econAct,h0RegimeOn,hxRegimeOn,numBootStep2);
switch_h0 = zeros(4,1);

for pVersion = 0:3
    if pVersion == 0
        switch_hx = [1 0 0 0;
            0 0 0 0;
            1 1 1 1;
            1 1 1 1];
    elseif pVersion == 1
        switch_hx = [1 1 1 1;
            1 1 1 1;
            0 0 0 0;
            0 0 0 0];
    elseif pVersion == 2
        switch_hx = [1 1 1 1;
            1 1 1 1;
            1 1 1 1;
            0 0 0 0];
    elseif pVersion == 3
        switch_hx = [1 1 1 1;
            1 1 1 1;
            0 0 0 0;
            1 1 1 1];
    elseif pVersion == 4
        switch_hx = [1 1 1 1;
            1 1 1 1;
            1 1 1 1;
            1 1 0 1];
    end
    [theta2Step3,AVarTheta2Step3] = step2SR_threshold_macro_spec(AVarTheta11Step3,outputStep3,bandwidthTstep2,setupStep3.econAct,switch_h0,switch_hx,numBootStep2);
    
    % Enforce consistency of sigma between Q and P
    posTheta12Step2 = zeros(1,length(namesTheta12));
    for i=1:length(namesTheta12)
        theta2Step3.(namesTheta12{i}) = theta2Step2.(namesTheta12{i});
        posTheta12Step2(1,i) = find(strcmp(AVarTheta2Step3.names,namesTheta12(i)));
    end
    AVarTheta2Step3.VarRobust(posTheta12Step2,posTheta12Step2) = AVarTheta2Step2.VarRobust(posTheta12Step2,posTheta12Step2);
    
    mprStep3 = getMPR(outputStep3,theta2Step3,AVarTheta2Step3,nm*2+nn,h0RegimeOn,hxRegimeOn);
    testRegimeStep3 = waldTestRegimeShift(theta2Step3,AVarTheta2Step3,2*nm+nn);
    
    %% Do the yield curve decomposition
    if h0RegimeOn == 0 && hxRegimeOn == 0
        yieldTP = yieldCurveDecom(setupStep3,outputStep3,theta2Step3);
    else
        numSimTP = 1D4;
        yieldTP  = yieldCurveDecom_threshold(setupStep3,outputStep3,theta2Step3,numSimTP);
    end
    
    % mean(yieldTP.termPremia(end,:))
    % std(yieldTP.termPremia(end,:))
    %% Simulating the model
    thresholdLevel          = AVarTheta2Step3.thresholdLevel;
    momentsData             = getMomentsData(loadData,outputStep3,thresholdLevel);
    momentsModel            = simulator_threshold(theta2Step3,outputStep3,thresholdLevel,numSim);
    %momentsModel_finite     = simulator_threshold_finite(theta2Step3,outputStep3,thresholdLevel,1000);
    expectedExcessReturn.h1 = getExpectedExcessReturns(loadData,setupStep3,outputStep3,theta2Step3,1);
    expectedExcessReturn.h3 = getExpectedExcessReturns(loadData,setupStep3,outputStep3,theta2Step3,3);
    
    %% Do plots
    plotFactors(setupStep3,outputStep3,AVarTheta11Step3)
    plotFit(setupStep3,outputStep3)
    plotUnconditionalMoments(momentsModel,momentsData)
    plotForecastRegressions(momentsModel,momentsData)
    plotYieldCurveDecom(yieldTP,setupStep3)
    
    %% Displaying the results
    outTheta1          = displayResults(theta1Step1,AVarTheta1Step1.VarRobust);
    outTheta11         = displayResults(theta11Step3,AVarTheta11Step3.VarRobust);
    theta2Step3Reduced = reduceTheta2(theta2Step3,AVarTheta2Step3.names);
    outTheta2Reduced   = displayResults(theta2Step3Reduced,AVarTheta2Step3.VarRobust);
    
    disp('optimum at step 1')
    printToScreen(struct2array(theta1Step1));
    disp('optimum at step 3')
    theta1Step3Start  = theta11Step3;
    for i = 1:size(namesTheta12,1)
        theta1Step3Start.(namesTheta12{i}) = setupStep3.calibrateTheta1.(namesTheta12{i});
    end
    printToScreen(struct2array(theta1Step3Start));
    
    %% Saving the results in a mat file
    if modelNum == 1
        modelName = ['QTSM_nn0_ReStudPhiQ_APversion',num2str(pVersion),'_'];
    elseif modelNum == 2
        modelName = ['QTSMr_nn0_ReStudPhiQ_APversion',num2str(pVersion),'_'];
    end
    filename = [modelName,num2str(startYear),'_SE',num2str(epsValue),'_h0RegimeOn',num2str(h0RegimeOn),'_','hxRegimeOn',num2str(hxRegimeOn),...
        '_numBootStep2_',num2str(numBootStep2)];
%     save([pwd,'\Results\',filename,'.mat'],'QStep1','theta1Step1','outputStep1','AVarTheta1Step1',...
%         'theta2Step2','AVarTheta2Step2',...
%         'theta11Step3','outputStep3','setupStep3','AVarTheta11Step3','theta2Step3','AVarTheta2Step3',...
%         'momentsModel','momentsData','mprStep3','testRegimeStep3','yieldTP','expectedExcessReturn',...
%         'switch_hx','switch_h0')
end
